General Relativistic Simulations of Black Hole-Neutron Star Mergers: 

Effects of black-hole spin 
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Black hole-neutron star (BHNS) binary mergers are candidate engines for generating both short- 
hard gamma-ray bursts (SGRBs) and detectable gravitational waves. Using our most recent confor- 
mal thin-sandwich BHNS initial data and our fully general relativistic hydrodynamics code, which 
is now AMR-capable, we are able to efficiently and accurately simulate these binaries from large 
separations through inspiral, merger, and ringdown. We evolve the metric using the BSSN formu- 
lation with the standard moving puncture gauge conditions and handle the hydrodynamics with a 
high-resolution shock-capturing scheme. We explore the effects of BH spin (aligned and anti-aligned 
with the orbital angular momentum) by evolving three sets of initial data with BH:NS mass ratio 
^ = 3: the data sets are nearly identical, except the BH spin is varied between a/MsH = —0.5 
(anti-aligned), 0.0, and 0.75. The number of orbits before merger increases with a/MeH, as ex- 
pected. We also study the nonspinning BH case in more detail, varying q between 1, 3, and 5. We 
calculate gravitational waveforms for the cases we simulate and compare them to binary black-hole 
waveforms. Only a small disk (< O.OIM©) forms for the anti-aligned spin case (a/MsH = —0.5) and 
for the most extreme mass ratio case (^ = 5). By contrast, a massive (Mdisk ~ O.2M0), hot disk 
forms in the rapidly spinning (a/MeH = 0.75) aligned BH case. Such a disk could drive a SGRB, 
possibly by, e.g., producing a copious flux of neutrino-antineutino pairs. 
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I. INTRODUCTION 

Mergers of black hole-neutron star (BHNS) binaries 
are expected to be among the most promising sources of 
gravitational waves detectable by ground-based laser in- 
terferometers like LIGO VIRGO SH, GEO 0, 
TAMA [6, 7], and AIGO [s], as well as by the proposed 
space- based interferometers LISA [9.] and DECIGO 0. 
Theoretical models indicate that a neutron star-neutron 
star (NSNS) 0, [TlEl^^ or black hole- 
neutron star (BHNS) pjl, [isl, M, llO, [IH, [22] binary 
merger may result in a hot, massive disk around a black 
hole, whose temperatures and densities could be suffi- 
cient to trigger a short-hard gamma-ray burst (SGRB). 
Indeed, SGRBs have been repeatedly associated with 
galaxies with extremely low star formation rates (see [23| 
and references therein for a review), indicating that the 
source is likely to involve an evolved population, rather 
than main sequence stars. The number of detectable 
BHNS mergers in the observable universe is still an open 
question. The uncertainties arise from some aspects of 
population synthesis calculations that are only partially 
understood. The estimated event rate of BHNS mergers 
observable by an Advanced LIGO detector typically falls 
in the range 7^ - 1 - 100 yr"^ 0. 

Motivated by the significance of BHNS binaries for 
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both GRB and gravitational- wave physics, many simula- 
tions of BHNS systems have been performed over the past 
years. In most dynamical simulations of BHNS binaries 
to date, the self-gravity of the neutron star (NS) and/or 
the tidal gravity of the black hole (BH) are treated in 
a Newtonian or post-Newtonian framework (see, e.g., 
m, m, 113, m, [il; see also [30], who provide fully rela- 
tivistic head-on collision simulations). 

Our first BHNS merger simulations were performed in 
an approximate relativistic framework [l^, |Q , which as- 
sumed conformal flatness of the spatial metric through- 
out the evolution (see |Q). In particular, we evolved 
extreme mass ratio {q = Mbh/^ns^I) initial data de- 
veloped earlier by our group [33, 34]. Though this simula- 
tion technique only allows for crude estimates, we found 
that mergers of irrotational BHNS binaries may lead to 
disks with masses up to 0.3 M©, with sufficient heating 
to emit the neutrino fluxes required to launch an SGRB. 

To date the only self-consistent, general relativistic, 
dynamical simulations of BHNS inspiral and coalescence 
have been those of Shibata et al [iO, Hi], [22|, |3^, Eti- 
enne et al. [36] (hereafter Paper I), and Duez et al. [37| . 
In many earlier calculations, especially those with initial 
mass ratios g ^ 3, significant disks are formed after the 
NS is disrupted, and, for very stiff nuclear equations of 
state (EOSs), the core of the NS may survive the initial 
mass transfer episode and remain bound. These findings 
contrast with some semi-analytic relativistic arguments, 
which suggest that disks with appreciable masses are dif- 
ficult to form via the merger of BHNS binaries [38] . 

In their first set of fully relativistic BHNS simulations, 
Shibata and Uryu found disk masses in the range of 0.1 
- 0.3 Mq for corotating neutron stars [lO, [2l|. Subse- 
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quently, Shibata and Taniguchi ([22|], hereafter ST) found 
smaller disk masses of 0.04 - 0.16 Mq for irrotational 
neutron stars, which are more realistic. However, our 
own fully general relativistic simulations of irrotational 
BHNSs in Paper I suggest that the disk mass is no more 
than 0.04 Mq. Such a small disk mass, if true for generic 
cases, might disfavor BHNS mergers as possible central 
engines for SGRBs. 

Recently, Yamamoto, Shibata, and Taniguchi pre- 
sented a new code which implements an adaptive mesh 
refinement (AMR) algorithm [35]. In one of their test 
simulations of BHNS binaries, they use the same initial 
data as in ST, but they find that the disk mass is less 
than 0.001 Mq. This result disagrees with that in ST, 
but is consistent with our result in Paper I. They at- 
tribute the discrepancy to different grid structures and 
computational methods. 

Most recently, Duez et al. developed a new fully gen- 
eral relativistic hydrodynamics code [37] that evolves the 
generalized harmonic formulation of Einstein's equations 
using a pseudospectral method, and the equations of hy- 
drodynamics via shock-capturing, finite-difference tech- 
niques. They use this code to perform a simulation of 
a single equal-mass BHNS and find that the disk mass 
is less than O.O3M0, which agrees with our result in Pa- 
per I. 

As we report in this paper, our code has been modified 
to use the moving-box Carpet [39] AMR infrastructure. 
With AMR, our code is now capable of performing BHNS 
simulations with resolutions equivalent to Paper I at only 
about ~ 1% the total computational cost. This enables 
us to simulate strong-field regions at higher resolution 
while simultaneously extending our grid's outer bound- 
aries much farther into the wavezone for more accurate 
wave extraction. We also evolve initial data at larger 
binary separations than in Paper I. Some of the BHNS 
initial data evolved in this paper are also new, as the BHs 
are now spinning. All of our initial data are constructed 
in the conformal thin-sandwich (CTS) formalism. The 
BH spin is incorporated by imposing boundary condi- 
tions on the BH horizon as in [40, 41]. The initial NS is 
irrotational and is modeled by an n = 1 polytropic EOS. 

In this paper, we present a new set of BHNS simu- 
lations that probes how the initial BH spin and the bi- 
nary mass ratio affect the dynamics and outcome of the 
merger. In particular, we study mergers of nonspinning 
BHNS binaries with mass ratios g = 1, 3, and 5. For 
the mass ratio (7 = 3, we study three cases in which the 
BH spin parameter (a = Jbh/^bh) —0.5 (counter- 
rotating), and 0.75. As in Paper I, we focus on binaries 
in which the NS tidally disrupts before reaching the in- 
nermost stable circular orbit (ISCO) and plunging into 
the BH. These systems are the most likely to yield a sig- 
nificant disk after merger. We thus do not evolve extreme 
mass-ratio binaries in this paper. 

Starting with approximately the same initial orbital 
angular velocity 1^, we find that when the BH is spinning, 
the merger time is delayed when the BH spin aligns with 



the orbital angular momentum, and is shortened when 
the BH spin is anti-aligned with the orbital angular mo- 
mentum. For example, evolutions of initial data with 
MVt ^ 0.033 [M is the initial Arnowitt-Deser-Misner 
(ADM) mass of the system], require about 4.25 orbits be- 
fore merger for a = —0.5, but 6.5 orbits before merger for 
a = 0.75. This finding is expected, since for large, aligned 
spins there is more angular momentum that must be ra- 
diated away to bring the binary to merger. This result 
is consistent with similar findings in the case of BHBHs 
(see e.g. jH). 

After the merger, we find a remnant disk with rest mass 
between < O.OIM© and ^ 0.2Mq (assuming the NS has a 
rest mass of Mq = I.4M0). When the initial BH is non- 
spinning, the disk mass is ~ O.O6M0 for = 3, ~ O.O3M0 
for g = 1 and tiny (^ O.OIM©) for q = 5. For a fixed 
mass ratio g = 3, the disk mass is ^ 0.2Mq for a = 0.75 
and tiny (< O.OIM©) for d = —0.5. The increase in disk 
mass with increasing BH spin agrees qualitatively with 
the semi-relativistic, smoothed particle hydrodynamics 
(SPH) simulations of large mass-ratio {q ^ 10) BHNS 
mergers reported in [29]. 

We find that the remnant disks are hot. A rough 
estimate suggests that the temperature in the disks is 
T - 5 X lO^^K. The results of BH disk simulations in [H 
suggest that the disk could generate a gamma-ray energy 
E ~ 10^^-lO^^erg from neutrino-antineutrino annihila- 
tion, which is promising for BHNS mergers as plausible 
central engines for SGRBs. 

We compute gravitational waveforms and the effective 
wave strain in the frequency domain for our mergers. As 
in Paper I, we find measurable differences between BHNS 
waveforms and those produced by BHBH mergers within 
the Advanced LIGO band. These differences appear at 
frequencies at which the NS is tidally disrupted and ac- 
creted by the black hole. Extracting the masses and spins 
of the objects from the inspiral data, the merger wave- 
forms could provide information about the radius of the 
NS, which in turn may be used to constrain the NS equa- 
tion of state. 

In contrast with Paper I, we find substantial disks in 
two of our nonspinning BH cases. The reason is that in 
Paper I, we imposed limits on pressure to prevent spu- 
rious heating/cooling in the low-density artificial "atmo- 
sphere". We find that imposing these artificial limits 
spuriously removes angular momentum, especially dur- 
ing the NS tidal disruption phase, thereby suppressing 
disk formation (see Section [VD] for further details). 

As we cautioned in Paper I, our findings are funda- 
mentally limited by the physics we do not model in these 
preliminary simulations (e.g., neutrino transport, mag- 
netic fields, etc), as well as by our choice of a simple 
F-law EOS to model the (cold and hot) nuclear EOS. 
Disk masses may depend sensitively on the initial struc- 
ture of the neutron star, especially the low-density outer 
regions, which is determined by the adopted EOS. So all 
BHNS merger results should be viewed in light of these 
caveats. Our main focus has been to establish that re- 
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liable simulations in full general relativity can now be 
launched to address these important issues. 

This paper is organized as follows. In Sees. HIl and lllli 
we briefly outline the basic equations and their specific 
implementation in our general relativistic, hydrodynam- 
ics scheme. Here we also provide an overview of our ini- 
tial data, gauge conditions, matter evolution equations, 
and diagnostics. Sec. IIVI presents code tests designed 
to validate our new AMR-based scheme. In Sec. [Vl we 
review the results of our BHNS merger simulations. In 
Sec. IVII we summarize our findings and comment on fu- 
ture directions. 



II. BASIC EQUATIONS 

The formulation and numerical scheme for our simu- 
lations are basically the same as those already reported 
in [H, to which the reader may refer for details. 
Here we introduce our notation, summarize our method, 
and point out the latest changes to our numerical tech- 
nique. Geometrized units {G = c = 1) are adopted, 
except where stated explicitly. Greek indices denote all 
four spacetime dimensions (0, 1, 2, and 3), and Latin 
indices imply spatial parts only (1,2, and 3). 

We use the 3+1 formulation of general relativity and 
decompose the metric into the following form: 

ds^ = -a^dt^ + -fij{dx' + P'dt){dx^ + P^dt) . (1) 

The fundamental variables for the metric evolution 
are the spatial three-metric jij and extrinsic curva- 
ture Kij. We adopt the Bau mgar te-Shapiro-Shibata- 
Nakamura (BSSN) formalism [45|,[4^ in which the evolu- 
tion variables are the conformal exponent (j) = ln(7)/12, 
the conformal 3-metric jij = e~^^^ij^ three auxiliary 
functions = — 7*-^ ,j, the trace of the extrinsic curva- 
ture and the tracefree part of the conformal extrinsic 
curvature Aij = e~^^{Kij — ^ijK /?>). Here, 7 = det(7ij). 
The full spacetime metric g^j, is related to the three- 
metric 7^^, by 7^^ = g^j, + n^n^,, where the future- 
directed, timelike unit vector normal to the time slice 
can be written in terms of the lapse a and shift (3^ as 
= 1,-/5*). The evolution equations of these 
BSSN variables are given by Eqs. (9)-(13) in Paper I. 

It has been suggested that evolving x — or 
W = e~'^^ instead of (j) gives more accurate results in 
BHBH simulations (see e.g. [13, S, |49| ) . We have tried 
all three techniques. We find that in the presence of hy- 
drodynamic matter, the evolution of the variable (j) yields 
a more stable evolution inside the apparent horizon (AH) 
near the puncture, leading to better rest-mass conserva- 
tion (see Sec. IIII B]) . We therefore adopt (/)-evolution in 
all of our BHNS simulations. 

It has been suggested that Kreiss-Oliger dissipation is 
sometimes useful to add in the BSSN evolution equations 
outside the AH to reduce high-frequency numerical noise 
associated with AMR refinement interfaces [50] . We have 



also tried this technique, and found that the dissipation 
results in slightly smaller Hamiltonian and momentum 
constraint violations. We typically do not use it in our 
BHNS simulations. However, in one case (Case E, see 
Sec. lYl), we found that the presence of hydrodynamic 
matter causes the conformal related metric 7^^ to lose 
positive definiteness near the puncture, which is unphys- 
ical. This behavior eventually causes the code to crash. 
We find that adding Kreiss-Oliger dissipation, coupled 
with a high but finite pressure ceiling deep inside the 
AH, solves this problem. 

We adopt standard puncture gauge conditions: an ad- 
vective "1+log" slicing condition for the lapse and a 
"Gamma- freezing" condition for the shift [5l|. Thus we 
have 

d^a = -2aK , (2) 
dop' = (3/4)5* , (3) 
doB' = dor - r]B' , (4) 

where do = dt — /3^dj. The 77 parameter is set to 
^ 0.55/M at the beginning of all simulations, where M is 
the initial ADM mass of the BHNS binary. If the value 
for 77 is fixed at this value throughout the simulation, 
both Hamiltonian and momentum constraint violation 
increase to marginally unacceptable levels (^ 10%, as 
measured by Eqs. (40) and (41) in Paper I) outside the 
BH during the merger phase, especially in the high BH 
spin cases. By simply doubling 77 to ~ 1.1 /M early in the 
inspiral phase (specifically, at t 55M) and fixing it at 
that value significantly reduces constraint violations out- 
side the BH during merger. However, late in the merger 
stage of Case B, we find that Hamiltonian constraint vi- 
olation increases rapidly (^ 10%) outside the BH, unless 
we enable Kreiss-Oliger dissipation as well. 

The fundamental matter variables are the rest-mass 
density po, specific internal energy e, pressure P, and 
four- velocity u^. We adopt a F-law EOS P = {T — l)po€ 
with r = 2, which reduces to an n = 1 poly tropic law for 
the initial (cold) neutron star matter. The stress-energy 
tensor is given by 

T^u = pohu^u^ + P^^^ , (5) 

where h = 1 -\- e -\- P/ po is the specific enthalpy. In our 
numerical implementation of the hydrodynamics equa- 
tions, we evolve the "conservative" variables p*, 5^, and 
f. They are defined as 

p* = -Vrpo^/i^^ , (6) 

f = ^T^,n^rf - p, . (8) 

The evolution equations for these variables are given by 
Eqs. (21)-(24) in Paper 1. 
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III. NUMERICAL METHODS 
A. Initial data 

Our initial data are constructed by solving Einstein's 
constraint equations in the GTS formalism, which allows 
us to impose an approximate helical Killing vector by set- 
ting the time derivatives of the conformally related metric 
jij to zero. We model the NS as an irrotational n = 1 
polytrope, and impose the black hole equilibrium bound- 
ary conditions of Cook and Pfeiffer [40] on the black hole 
horizon. Details of this method can be found in [52]. Our 
initial data here differ from those described in [52] only 
in terms of the black hole spin, which can be specified 
by virtue of a free (vectorial) parameter ^7! that appears 
in the boundary conditions of [40]. In \E2\ we adopted 
the method in [41] by iterating over this parameter until 
the quasilocal spin Jbh of the black hole vanishes, mak- 
ing both binary components irrotational. Here we also 
consider rotating black hole configurations, for which we 
iterate over l^J: until the black hole spin equals certain 
specified values. We focus on black hole spins that are 
aligned or anti-aligned with the orbital angular momen- 
tum, and consider a = Jbh/^bh ~ ^-^^ (aligned), 0.0 
(nonspinning), and —0.5 (anti-aligned). 

The initial data are calculated using the Lorene spec- 
tral methods numerical libraries \E§\ . To map these spec- 
tral configurations onto our non-spectral simulation grid, 
we first construct our numerical grid and record the po- 
sitions of each point in physical coordinates. Then we 
evaluate the field and hydrodynamic quantities at these 
points based on their spectral coefficients [54]. Finally, 
the excised BH region is filled with constraint-violating 
initial data, using the "smooth junk" technique we devel- 
oped and validated in [55] (see also [5^,(53]). In partic- 
ular, we extrapolate all initial data quantities from the 
BH exterior into the interior with a 7th order polynomial, 
using a uniform stencil spacing of Ar ^ 0.3rAH- 

All of the NSs considered in this paper have a com- 
paction of C = Mns/^ns = 0.145, where Mns is the 
ADM mass and R^s is the (circumferential) radius of the 
NS in isolation. Since we model the NS with an n = 1 
(r = 2) polytropic EOS, the rest mass of the NS, Mq, 
scales with the polytropic constant k, as Mq oc z^^/^. For 
a NS with compaction C = 0.145, we find the ADM mass 
for the isolated NS to be Mns = 1.30Mq{Mo/IAMq), 
with an isotropic radius i^iso = 11.2km(Mo/1.4M0) 
and circumferential (Schwarzschild) radius of R^s = 
13.2km(Mo/1.4M0). The maximum rest-mass density 
of this NS is /)o,max = 9 X lO^'^g cm-3(1.4M0/Mo)^ 

B. Evolution of the metric and hydrodynamics 

We evolve the BSSN equations with fourth-order ac- 
curate, centered finite-differencing stencils, except on 
shift advection terms, where we use fourth-order accurate 
upwind stencils. We apply Sommerfeld outgoing wave 



boundary conditions to all BSSN fields, as in Paper I. Our 
code is embedded in the Cactus parallelization frame- 
work [5^, and our fourth-order Runge-Kutta timestep- 
ping is managed by the MoL (Method of Lines) thorn, 
with a Courant-Friedrichs-Lewy (CFL) factor set to 0.45 
in all BHNS simulations. We use the Carpet [39] infras- 
tructure to implement the moving-box adaptive mesh re- 
finement. In all AMR simulations presented here, we use 
second-order temporal prolongation, coupled with fifth- 
order spatial prolongation. The apparent horizon (AH) 
of the black hole is computed with the AHFinderDirect 
Cactus thorn [59|. 



We write the general relativistic hydrodynamics equa- 
tions in conservative form. They are evolved by a 
high-resolution shock-capturing (HRSC) technique Q 
that employs the monotonized central (MC) reconstruc- 
tion scheme [60j coupled to the Harten, Lax, and van 
Leer (HLL) approximate Riemann solver [61]. The 
adopted hydrodynamic scheme is second-order accurate 
for smooth flows, and first-order accurate when disconti- 
nuities (e.g. shocks) arise. To stabilize our hydrodynamic 
scheme in regions where there is no matter, we maintain 
a tenuous atmosphere on our grid, with a density floor 
Patm set equal to 10 times the initial maximum density 
on our grid. The initial atmospheric pressure Patm is set 
equal to the cold polytropic value Patm = ^Patm? where 
is the polytropic constant at t = 0. Throughout the evo- 
lution, we impose limits on the atmospheric pressure to 
prevent spurious heating and negative values of the inter- 
nal energy e. Specifically, we require Pmin < P < Pmax, 
where Pmax = lOhzp^ and Pmin = ^Po/^- Whenever P 
exceeds Pmax or drops below Pmin, we reset P to Pmax 
or Pmin, respectively. Applying these limits everywhere 
on our grid would artificially sap the angular momentum 
in the tidally disrupted NS, allowing matter to fall spuri- 
ously into the BH and thereby suppressing disk formation 
(see Sec.|V|). To effectively eliminate this spurious angu- 
lar momentum loss, we impose these pressure limits only 
in regions where the rest-mass density remains very low 
{po < lOOpatm) or deep inside the AH, where e^^ > V^^ax 
(see Sec. IVDI for more details). 



At each timestep, we need to recover the "primitive 
variables" po, P, and v'^ from the "conservative" variables 
p*, f, and Si. We perform the inversion as specified in 
Eqs. (57)-(62) of [4^, but with a slightly modified ana- 
lytic quartic solver from the GNU Scientific Library that 
outputs only the real roots. We use the same technique 
as in Paper I to ensure that the values of Si and f yield 
physically valid primitive variables, except we reset f to 
^o,max (where To,max is the maximum value of r ini- 
tially) when either Si or f is unphysical [i.e., violate one 
of the inequalities (34) or (35) in Paper I]. The restric- 
tions usually apply only to the region near the puncture 
or in the low-density atmosphere. 
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C. Diagnostics 

During the evolution, we monitor the Hamiltonian and 
momentum constraints calculated by Eqs. (40)-(43) of 
Paper I. Gravitational waves are extracted using both 
the Regge- Wheeler- Zerilli-Moncrief formalism and the 
Newman-Penrose Weyl scalar t/j^. The extraction tech- 
nique is summarized in Sec. IIIF of Paper 1. 

We also monitor the ADM mass and angular momen- 
tum of the system, which can be calculated by surface in- 
tegrals at a large distance (Eqs. (37) and (39) of Paper I). 
However, with our AMR grid, the resolution is rather low 
in regions very far from the binary, which causes our an- 
gular momentum measurement to suffer from numerical 
noise. This problem can be overcome by using Gauss's 
law to split the distant surface integral into a sum of in- 
tegrals over an inner surface and the volume between the 
inner and distant surfaces. The reasons are twofold: (1) 
the inner surface integrals are computed more accurately, 
as the inner surface resides in the region where numeri- 
cal resolution is higher, and (2) most of the contribution 
in the volume integral is from the strong-field domain, 
which is also in the inner (high-resolution) region. The 
equations we use to calculate ADM mass and angular 
momentum with minimal numerical noise are as follows 
(see e.g. Appendix A in [62] for a derivation): 

M = [ d^x (^''p+ -^4>'A,,A'^ - -^r^'^Fjik (9) 

Ji = ^^ij'^ j^<i'x[i,\A^n + \x^dnK (10) 
-^X^Akmdnfn + ^T^xiSn] 

Here 5 is a surface surrounding the BH, V is the vol- 
ume between the inner surface S and a distant surface, 

= e^, p = n^n^.T^'', Si = -n^ji^T^'', R is the Ricci 
scalar associated with 7^^, and Tijk are Christoffel sym- 
bols associated with jij . 

When hydrodynamic matter is evolved on a fixed uni- 
form grid, our hydrodynamic scheme guarantees that the 
rest mass Mq is conserved to machine roundoff error. 
This strict conservation is no longer maintained in an 
AMR grid, where spatial and temporal prolongation is 
performed at the refinement boundaries. Hence we also 
monitor the rest mass 

^0 = J P^d^x (11) 

during the evolution. Rest-mass conservation is also vio- 
lated whenever po is reset to the atmosphere value. This 



usually happens only in the very low-density atmosphere 
or deep inside the AH where accuracy is not maintained. 
The low-density regions do not affect rest-mass conser- 
vation significantly. However, this is not the case near 
the puncture, where the rest-mass density can be very 
high, especially after the NS is tidally disrupted and mat- 
ter falls into the BH and accumulates near the puncture. 
The evolution near the puncture is sensitive to the choice 
of conformal variable used in the BSSN evolution, as dis- 
cussed in Sec. [Hi We find that the (j) evolution leads to 
better rest-mass conservation than x or evolutions. 

We measure the thermal energy generated by shocks 
via the quantity K = P/Pcoid, where Pcoid = i^Po is the 
pressure associated with the cold EOS that characterizes 
the initial matter. The specific internal energy can be 
decomposed into a cold part and a thermal part: e = 
ecoid + eth with 

ecoid = - j Pco\dd{l/po) = f^Po"' • (12) 

Hence the relationship between K and eth is 

_ 1 P ^ r-i 

^th — e - ecoid — 7 tPo 

i — 1 Po i — 1 

= {K- l)ecoid . (13) 

For shock-heated gas, we always have K > 1 (see Ap- 
pendix [B]). 

IV. AMR CODE TESTS 

Our general relativistic magnetohydrodynamic 
(GRMHD) code has been thoroughly tested by passing 
a robust suite of tests. These tests include maintaining 
stable rotating stars in stationary equilibrium, reproduc- 
ing the exact Oppenheimer- Snyder solution for collapse 
to a black hole, and reproducing analytic solutions 
for MHD shocks, nonlinear MHD wave propagation, 
magnetized Bondi accretion, and MHD waves induced 
by linear gravitational waves 0, [H] . Our code has also 
been compared with Shibata & Sekiguchi's GRMHD 
code by performing simulations of identical mag- 
netized hypermassive neutron stars [H, [66| and of the 
magnetorotational collapse of identical stellar cores [l^ . 
We obtain good agreement between these two indepen- 
dent codes. Our code has also been used to simulate 
the collapse of very massive, magnetized, rotating stars 
to black holes [67]; merging BHBH binaries [55], BHNS 
binaries (Paper I), magnetized NSNS binaries [68,]; and 
relativistic hydrodynamic matter in the presence of 
puncture black holes [69^. Recently, our code has been 
generalized to incorporate (optically thick) radiation 
transport and its feedback on MHD fluids [7Q| . 

All of the above-mentioned tests and simulations were 
performed on grids with uniform spacing. In some of 
the simulations, we utilized the multiple-transition fish- 
eye transformation [71] so that a uniform computational 
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grid spacing corresponds to physical coordinates with 
spatially varying resolution. Recently, we have modi- 
fied our code so that we can use the moving-box AMR 
infrastructure provided by Carpet [39]. To test our new 
code, we have performed shock-tube tests and 3+1 simu- 
lations of linear gravitational waves, single stationary and 
boosted puncture BHs, puncture BHBH binaries, stable, 
rapidly and differentially rotating relativistic stars, and 
relativistic Bondi accretion onto a Schwarzschild BH. All 
of our 3+1 AMR code tests were performed assuming 
equatorial symmetry (i.e., symmetry about the z = 
orbital plane), which we assume in all BHNS evolutions 
presented in this paper. In many of the tests we have an- 
alytic solutions with which to compare. For the BHBH 
tests, we compare our results with the numerical results 
reported in the literature. We have confirmed that our 
new code can evolve these systems reliably with AMR. 
Below we briefly report results from both our BHBH and 
differentially rotating relativistic star AMR code tests. 



A. Binary black holes 

In this subsection we summarize an important testbed 
for calibrating the vacuum sector of our code: BHBH 
merger simulations. We provide details on the setup of 
our initial data and numerical grid, as well as on the 
adopted numerical techniques and simulation results, in 
the tables and figures of Appendix [Al Here we point out 
some highlights. 

We first simulate an equal-mass BHBH system. We 
use puncture initial data [72, [t^, [zl] with initial binary 
separation D = 9.89M, where M is the ADM mass of 
the binary. The initial "bare" mass and momentum of 
each puncture are chosen according to [75] to make the 
orbit quasicircular (see also [76]). This initial configura- 
tion is close to the R4 run presented in (HO, [ij • We use 
8 refinement levels centered at each BH and place the 
outer boundary at 323. 4M. The resolution in the inner- 
most refinement level is M/32, so that there are about 
15 grid points across the diameter of each black hole's 
AH initially. This number rapidly increases during the 
early inspiral to a fixed value of about 24 grid points 
due to our gauge choice, and increases again to about 
40 grid points, post-merger. We find that the binary in- 
spiral orbit, gravitational waveform, energy and angular 
momentum emitted by gravitational waves, and the fi- 
nal black hole's spin parameter are very similar to those 
reported by others in [10, [ij • 

To further validate our numerical results, we compute 
the quantities 

SE = {M - Mf - AEgw)/M , (14) 
SJ = {J-Jf-AJGw)/J , (15) 

where J is the ADM angular momentum of the initial 
binary, Mf and Jf are the ADM mass and angular mo- 
mentum of the final (merged) system, and A£^gw and 



A Jgw are the energy and angular momentum carried off 
by gravitational waves. We find that both Mf and Jf 
are very close to Mbh and Jbh, the mass and angular 
momentum of the final BH. We compute Mbh and Jbh 
from the irreducible mass Mirr and the ratio of proper 
polar to equatorial circumference Cr of the final black 
hole's horizon. Specifically, we first use Eq. (5.3) of [tI] 
to solve for the dimensionless spin d = Jbh/M^j^: 



where E{x) is the complete elliptic integral of the 
second kind. We next compute Mbh by Mbh = 

(Mirr/a)^2(l- Vl-a2) and Jbh by Jbh = clM^^. 
These formulae are derived for a Kerr spacetime. They 
are applicable to the merged BH as the spacetime ap- 
proaches Kerr at late times. Gravitational waves are ex- 
tracted via the Newman-Penrose scalar ?/^4. The quan- 
tities A£^GW and AJqw are computed by integrating 
Eqs. (50) and (52) of Paper I. Conservation of energy 
and angular momentum demands that SE = = SJ. 
Hence these parameters are good indicators of the accu- 
racy of the simulations, whenever a significant amount of 
energy and angular momentum are radiated during the 
evolution. We found SE ^ 4 x IQ-"^ and SJ^4x 10"^ 
in our equal-mass BHBH simulation. 

We next perform an unequal-mass BHBH simulation. 
The initial data are generated by the TwoPunctures 
code [79] with initial binary (coordinate) separation D ^ 
7M. The approximate initial momentum of each BH for 
a quasicircular configuration is calculated by the 3PN ex- 
pression using Eq. (65) of [47]. The "bare" masses of the 
punctures are adjusted so that the irreducible mass ratio 
of the BHs equals 3. The initial configuration approxi- 
mates one of the cases studied in [80|. We use 9 and 10 
refinement levels centered at the larger and smaller BH, 
respectively, with an outer boundary at 413. 9M. The 
resolution around the large and small BH is M/50 and 
M/lOO, respectively, so that there are initially about 36 
and 22 grid points across the diameter of the larger and 
smaller black holes' AH, respectively. We find that our 
computed values of A£^gw and AJqw agree with Ta- 
ble I of [80] to within 4%, and that the kick velocity 
agrees with Fig. 2 of [8l|. We note that exact agree- 
ment is not expected as our initial configurations are 
slightly different. Finally, we find excellent energy and 
angular momentum conservation: SE ^ —2 x 10~^ and 
SJ ^9x 10-^. 

These tests demonstrate our code's ability to evolve 
dynamical vacuum spacetimes containing moving BHs. 
We next study our code's ability to handle hydrodynam- 
ics in a strong- field spacetime with moving AMR boxes. 
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FIG. 1: The moving AMR grid in the rotating relativist ic star test. Shown here are the snapshots of the two innermost 
refinement boxes (the nested squares) and part of the third refinement level in the equatorial plane. The circle denotes the 
boundary of the star. 



B. Rotating relativist ic star 

In our BHNS simulations, the NS is placed inside a 
high resolution refinement box initially. However, when 
it is tidally disrupted by the black hole, the NS matter 
spreads into other refinement boxes, with different reso- 
lutions. The purpose of this test is to demonstrate the 
ability of our AMR code to maintain the equilibrium of 
a stable, differentially rotating relativistic star crossing 
refinement boundaries during the evolution, and verify 
second-order convergence as the resolution is improved. 

The star considered in this test is the same hypermas- 
sive neutron star studied in [82] , which is modeled by an 
n = 1 polytropic EOS. The mass of the star is 1.7 times 
greater than the maximum mass of a nonrotating neutron 
star with the same EOS, but has been shown to be dy- 
namically stable in previous simulations 0, [82|, [^] . The 
star is highly flattened because of its rapid rotation. It 
has an equatorial (coordinate) radius Re ^ 4.5M, and a 
polar radius Rp ^ 1.2M. The rapid rotation also causes 
a torus-like density distribution in which the maximum 
density occurs off-center. The central rotation period of 
the star is Prot = 38.4M. 

Our AMR grid contains four refinement levels that 



move in time. The side-lengths of the refinement boxes 
in the equatorial plane are 4.29M, 8.69M, 22. 9M, and 
45. 8M. Since we impose equatorial symmetry, the 
lengths in the 2:-direction are half of the sizes stated 
above. The coarsest grid is fixed and the outer boundary 
is placed at x, 7/ = ±46M and z = 43M. We initially cen- 
ter the star at the origin, and move the center of the four 
refinement boxes according to the following prescription: 

Xc = Ax sincot , He = Ay{l — cos cot) , = , (17) 

where A^^ Ay and uo are constants. We set the amplitudes 
Ax = 1.79M, Ay = 1.43M and the angular frequency 
of the moving center uj = 0.168/M (27r/cj = 0.98Prot). 
With these parameters, the bulk of the star crosses the 
innermost refinement boxes during the evolution (see 
Fig. [1]). The grid spacing doubles at each successive re- 
finement level. We evolve the star with three different 
resolutions, which we label low, medium and high resolu- 
tion runs. The grid spacing in the innermost refinement 
level is M/11.2, M/16.7, and M/22.2 for the low, medium 
and high resolution runs, respectively. 

Figure [2] shows the evolution of the rest mass Mq, ADM 
mass M, and angular momentum J of the star, calculated 
by Eqs. ([9|)-(pT]). Since there is no black hole in the space- 



FIG. 2: Evolution of rest mass Mo, ADM mass M, and an- 
gular momentum J of a rapidly rotating star using a moving 
AMR grid with low (black solid lines), medium (red dotted 
lines) and high (blue dashed lines) resolution. 

time, we calculate only the volume integrals of M and J 
over the entire spatial grid. As mentioned in Sec. IIII C[ 
strict rest-mass conservation is not expected when using 
an AMR grid. However, the deviations in Mq, M and 
J converge to zero at slightly higher than second order 
with increasing resolution. In addition, the amplitudes 
of the gravitational radiation decrease to zero at second 
order with increasing resolution. Gravitational radiation 
does not contribute significantly to the nonconservation 
of M and J observed in Fig. [2 AEqw/M < 10"^ and 
^Jgw/J ^ 10"^ for all three runs. 

Figure [3] demonstrates that the equilibrium profile of 
the rest-mass density is well-maintained in the high res- 
olution simulation, apart from a slight drift of the star's 
center. 



V. RESULTS 

A. Overview 

We perform a number of BHNS simulations with vary- 
ing initial configurations and numerical resolutions. Ta- 
ble [J provides an overview of our chosen initial configura- 
tions, and Table HIl specifies the AMR grid structure used 
in each case. The first three cases in Table |I] (Cases A, A- 
MSep and A-SSep) correspond to the same BHNS binary 
system tracked from different initial separations along a 
quasiequilibrium sequence. The objective of these sim- 
ulations is to demonstrate that the fate of the binaries. 



FIG. 3: Density profiles in the equatorial plane along the 
X-axis (upper graph) and ?/-axis (lower graph) for the high 
resolution run. The rest-mass density po is normalized to 
the initial maximum density, /?o,max- Solid (black), dotted 
(red) and dashed (blue) lines show the profiles at times t = 0, 
3.43Prot and 6.86Prot, respectively. 



TABLE I: Initial BHNS configurations. Here, the mass ratio 
is defined as ^ = Mbu/Mns , where Mbh and Mns are the 
ADM masses of the black hole and neutron star at infinite 
separation, M is the total ADM mass of the binary system, 
Do the initial binary coordinate separation, d = Jbh/M^^l 
is the spin parameter of the black hole (always aligned or 
anti-aligned with the orbital angular momentum), and Q is 
the initial orbital angular velocity. All neutron stars have the 
same nondimensional rest (baryon) mass Mo = Mq/k,^^'^ — 
0.15, which is ~ 83% the maximum rest mass of a nonrotating 
NS with the same n = 1 polytropic EOS. 



Case 


d 


q. 


Dq/M J/M'^ 




A 


0.00 


3.0 


8.81 


0.702 


0.0333 


A-MSep 


0.00 


3.0 


7.17 


0.668 


0.0441 


A-SSep 


0.00 


3.0 


5.41 


0.638 


0.0623 


B 


0.75 


3.0 


8.81 


1.096 


0.0328 


B-SSep 


0.75 


3.0 


5.53 


1.011 


0.0594 


C 


-0.50 


3.0 


8.81 


0.443 


0.0338 


D 


0.00 


5.0 


8.79 


0.518 


0.0333 


E 


0.00 


1.0 


8.61 


0.938 


0.0347 



and their merger waveforms, are insensitive to the ini- 
tial binary separation at which we begin the simulation, 
provided it is sufficiently large. For Case B-SSep, we 
evolve the system with three different resolutions (see 
Table [H]) to demonstrate convergence. Fixing the mass 
ratio dX q = ?> and initial orbital angular frequency at 
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TABLE II: Grid configurations used in our simulations. Here, A^ah denotes the number of grid points covering tfie diameter 
of the (spherical) apparent horizon initially, and A^ns denotes the number of grid points covering the smallest diameter of the 
neutron star initially. 



Case 


MO 


vjriiQ nieiaiCiiy yiiv uniLb oi ivi j 


IVIax. resolution 


A^AH A'ns 


A 


0.0333 


(196.7, 98.35, 49.18, 24.59, 12.29, 6.147, 3.073, 1.414 [1.660]) 


M/32.5 


41 


85 


A-MSep 


0.0441 


(197.0, 98.49, 49.24, 24.62, 12.31, 6.156, 3.078, 1.416 [1.662]) 


M/32.5 


40 


82 


A-SSep 


0.0623 


(197.3, 98.67, 49.33, 24.67, 12.33, 6.167, 3.083, 1.418 [1.665]) 


M/32.4 


40 


76 


B 


0.0328 


(210.2, 92.49, 46.24, 23.12, 11.56, 5.780, 2.890, 1.445 [1.642], 0.7554 [N/A]) 


M/60.9 


56 


80 


B-SSep-HRes 


0.0594 


(197.6, 98.79, 49.40, 24.70, 12.35, 6.174, 3.087, 1.544, 0.7718 [N/A]) 


M/64.8 


59 


78 


B-SSep-MRes 


0.0594 


(197.6, 98.79, 49.40, 24.70, 12.35, 6.174, 3.087, 1.544, 0.7718 [N/A]) 


M/47.9 


43 


57 


B-SSep-LRes 


0.0594 


(197.6, 98.79, 49.40, 24.70, 12.35, 6.174, 3.087, 1.544, 0.7718 [N/A]) 


M/41.5 


37 


50 


C 


0.0338 


(196.6, 98.31, 49.16, 24.58, 12.29, 6.145, 3.072, 1.413 [1.659]) 


M/32.5 


37 


85 


D 


0.0333 


(196.3, 98.12, 49.06, 24.53, 12.26, 4.415, 2.208, 1.104) 


M/48.9 


69 


84 


E 


0.0347 


(217.6, 108.8, 54.40, 27.20, 13.60, 6.801, 3.400, 1.496 [N/A], 0.7349 [N/A]) 


M/58.8 


47 


78 



There are two sets of nested refinement boxes: one centered on the NS and one on the BH. This column specifies the half 
side length of the refinement boxes centered on both the BH and NS. When the side length around the NS is different, we 
specify the NS half side length in square brackets. If there is no corresponding NS refinement box (as is the case when the NS 
is significantly larger than the BH), we write [N/A] for that box. 



TABLE HI: BHNS simulation results. Here Mdisk is the rest mass of the material outside the AH at the end of the simulation, 
a/ is the value of d computed by solving Eq. (|16|) at late times. The total energy and angular momentum carried off by the 
gravitational radiation are given by AEgw and A Jew, respectively. Vkick is the kick velocity due to recoil. Aorbits specifies 
the number of orbits in the inspiral phase before merger. 6E and 6 J measure the violation of energy and angular momentum 
conservation, as defined in Eqs. (|14|) and (|15|) . respectively. 



Case 


Mdisk/Mo 




AEgw/M AJgw/J 


V]^ick (km/s) 


A^orbits 


6E 


6J 


A 


3.9% 


0.559 


0.93% 


17.4% 


33 


4.5 


-0.02% 


2.2% 


A-MSep 


3.8% 


0.560 


0.81% 


13.2% 


33 


2.5 


-0.02% 


1.9% 


A-SSep 


2.9% 


0.557 


0.45% 


11.4% 


33 


1.75 


-0.02% 


-0.36% 


B 


15% 


0.881 


0.92% 


13.3% 


22 


6.5 


-0.04% 


2.8% 


B-SSep 


14.3% 


0.882 


0.60% 


6.4% 


79 


2.0 


-0.02% 


2.6% 


C 


0.8% 


0.331 


0.98% 


24.4% 


49 


3.25 


-0.01% 


1.3% 


D 


0.8% 


0.418 


0.98% 


19.2% 


73 


6.25 


-0.02% 


1.9% 


E 


2.3% 


0.851 


0.35% 


7.2% 


17 


2.25 


-0.01% 


0.9% 



MQ ^ 0.033, we then study the effects of black hole spin 
by choosing the spin parameter d between —0.5 and 0.75 
(Cases C, A and B). In all cases, the black hole spin 
is either aligned or anti-aligned with the orbital angular 
momentum. Finally, we study the dependence on the 
mass ratio by varying q between 1 and 5 (Cases E, A and 
D). 

In all simulations, AMR refinement levels are initially 
centered on the BH and NS. After t = 0, the levels cen- 
tered on the BH track the AH centroid, and the levels 
centered on the NS track the matter centroid defined by 

Xl^k^, (18) 

where V is the volume outside the innermost refinement 
box surrounding the BH. Figured shows snapshots of the 



inner refinement boxes during the Case A simulation. 

In all of our BHNS simulations, we find that the inte- 
grated Hamiltonian constraint violation [as measured by 
Eq. (40) in Paper I] is 0.1%-0.6% outside the BH during 
the inspiral phase. During merger, it jumps to a peak 
value of l%-2% in Cases A, C, D, and E and to 8.8% 
in Case B. Post-merger, the violation stays roughly con- 
stant (Case E) or gradually decreases (all other cases) 
until the simulation is stopped. The integrated momen- 
tum constraint violation [as measured by Eq. (41) in Pa- 
per I] outside the BH generally oscillates around 1%- 
2%, except for a few peaks of 3%-6% during inspiral. 
We measure energy and momentum conservation by the 
quantities SE and SJ defined in Eqs. ([Hj) and (p!5]) . We 
find that dE ^ 0.02% and 6 J = 1-2% (see Table [TTll) in 
all cases. 

We also monitor the rest-mass conservation by com- 
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FIG. 4: AMR grid setup during a typical BHNS simulation, as observed in the equatorial plane. Only the two innermost 
refinement boxes are shown for clarity. The black region denotes the apparent horizon of the black hole. Neutron star density 
contours are drawn according to po = po,maxlO~°"^^^~° °^ 0=0, l,--- 12)- The maximum rest-mass density of the initial NS is 
/^po,max = 0.126, or po,max = 9 X 10^"^ {1 AM Q / Mof g cm-\ 



puting the quantity 

AMo = [Mo - Mo(0)]/Mo(0) , (19) 

where Mo(0) is the NS rest mass computed at t = 0. 
Rest-mass conservation requires AMq = at all times. 
Most of the violation occurs inside the BH apparent hori- 
zon when the NS matter falls into the BH. We find that 



AMo ^ 0.06-0.08% during the inspiral phase (for all 
cases except E, which has AMq ^ 0.002%). During and 
after the merger phase, we find large "glitches" in Mq 
inside the AH due to inaccuracies in that region, which 
causes large AMq (^3%), but these are easily accounted 
for and neglected. We find that AMq is 0.1%-0.15% af- 
ter the inspiral phase. This result is important because 
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FIG. 5: The BH's irreducible mass Mirr and Cp/Ce versus 
time, for Case A. 

in some cases we find a low-mass (~ O.OIMq) disk sur- 
rounding the BH. These results can only be trusted if 
the rest-mass violation outside the apparent horizon due 
to numerical error is much smaller than the disk mass, 
which is true in all of our simulations. 

Figure [5] shows the BH's irreducible mass Mirr and 
Cr = Cp/Ce as a function of time for Case A, where Cp 
and Ce are the proper polar and equatorial circumfer- 
ences, respectively. Notice that these two quantities re- 
main almost constant at early times until the NS matter 
falls into the BH. They approach constant values again 
at late times as the system settles into a quasistation- 
ary state. We find the same qualitative behavior for all 
of our BHNS simulations. Given Cr, we define a us- 
ing Eq. (p!6|) . When the BH is well- approximated by the 
Kerr metric (e.g., the early inspiral phase, at large sep- 
arations), a equals the BH's spin parameter Jbh/M|jj. 
In its final stage, the system consists of a BH surrounded 
by a disk. The spacetime near the BH therefore devi- 
ates from Kerr, and a does not need to coincide with the 
value of the BH spin parameter defined in, e.g. the iso- 
lated horizon formulation [84]. Nevertheless, a is still a 
reasonable measure of the BH's spin when the disk is not 
very massive compared to the BH's mass. Table Hill lists 
the final value of a in each of our BHNS simulations. 



B. Convergence studies 

To demonstrate that the resolution used in our sim- 
ulations is sufficiently high, we perform a convergence 
test. It is known that as the BH spin is increased, higher 



FIG. 6: Numerical convergence of gravitational wave signals 
for three B-SSep simulations. The highest grid resolutions 
are M/41.5, M/47.9, and M/64.8, so that initially the AH 
diameter is covered by 38, 44, and 59 gridpoints, and the NS 
is covered by 58, 67, and 91 gridpoints, respectively. In the 
top panel, we show the real part of the / = m = 2 component 
of -04 for the three waveforms. In the bottom panel, pairwise 
differences between the waveforms are plotted and rescaled to 
demonstrate second-order convergence. 



resolution is required for accurate BHBH evolutions (see 
e.g. [85]). Hence we perform a convergence test for the 
highest spinning case with a = 0.75. We use three differ- 
ent grid resolutions with close initial binary separation 
(Cases B-SSep-HRes, B-SSep-MRes and B-SSep-LRes in 
Table ini), which is adequate for the purpose of a conver- 
gence test. The top panel in Fig. [6] shows the real part of 
?/^4^, the I = m = 2 mode of ?/^4, computed with three dif- 
ferent resolutions, showing good agreement. The bottom 
panel in Fig. [6] compares the pairwise differences between 
the waveforms and we see that the waveforms converge 
at second-order accuracy. 



C. Effects of different initial binary separations 

In addition to the convergence test, we need to ver- 
ify that our initial binary separation is sufficiently large. 
That is, starting the simulation from any larger initial 
separation should give the same outcome after a simple 
time (and, for gravitational waves, a possible phase) shift. 

Figure [71 shows the rest mass of NS matter outside the 
AH as a function of time for Cases A, A-MSep, and A- 
SSep. These cases correspond to initial binaries along the 
same quasiequilibrium sequence but with different initial 
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FIG. 7: Rest mass of NS matter outside the BH versus time 
for different initial separations for the q — 3, nonspinning 
BH case (Cases A, A-MSep, and A-SSep). Here, the time 
coordinate is shifted by t25, where tAr is the time at which 
N% of the NS rest mass has fallen into the apparent horizon. 



separations. The time origin is shifted by ^25, where In 
denotes the time at which N% of the NS rest mass has 
fallen into the AH. We see that the results in Cases A and 
A-MSep agree very well, but the result for A-SSep devi- 
ates from the other two. Figure [8] shows the gravitational 
waveforms Re('04^) emitted in Cases A and A-MSep. We 
again see that the waveforms agree very well after a sim- 
ple time and phase shift. These results indicate that ini- 
tial separations corresponding to those in Cases A-MSep 
and A are sufficiently large. For almost all other simu- 
lations, the initial value of [M/{q + 1)]1] is similar to or 
smaller than that for Case A {[M/{q + l)]n < 0.0083), 
suggesting that the initial binary separation should be 
sufficiently large in all of these cases. The only excep- 
tion is Case E, for which [M/{q + ^ 0.017 initially. 
This value should be compared with that of Case A-SSep, 
which had [M/{q-\- 1)]^ ~ 0.016 initially. Therefore, our 
results for Case E may be inaccurate, since the initial 
separation may be too small. From an astrophysical per- 
spective, our Case E may also be the least relevant, since 
the formation of a BH with a NS of equal mass seems 
very unlikely. 



D. Effects of pressure ceiling 

Figure [3 shows that a disk of rest mass 0.04Mo re- 
mains outside the AH at the end of simulation, where 
Mo is the rest mass of the initial NS. This finding dis- 



FIG. 8: Comparison of gravitational waveforms for Case A 
(black solid line) and A-MSep (blue dashed line). 



agrees with the result reported in Paper I. The reason 
for this discrepancy is due to artificially imposed limits 
on pressure P that were applied everywhere in the fluid 
in Paper I. Here, we apply these limits only in regions 
where the density is tiny (i.e. inside the artificial "atmo- 
sphere") or where the conformal factor tjj is large (i.e., 
deep inside the AH), as discussed in Sec. IIIIBl Setting 
the pressure ceiling everywhere has little effect on the 
inspiral dynamics. However, when the NS is tidally dis- 
rupted and forms a relatively low-density "tail", shock 
heating becomes a significant effect, causing P to exceed 
^max (see Appendix [B]) . In our original formulation, we 
reset P to Pmax- Artificially lowering the pressure in this 
way forces the low-density material to lose a significant 
amount of torque and angular momentum and fall into 
the BH. Even though the density of the tail is low, it is 
still many orders of magnitude larger than the artificial 
atmospheric density, and hence has substantial influence 
on the dynamics of the merging NS. 

To more accurately model the physics of the merger 
phase and maintain a stable evolution inside the BH, 
we now impose the limits on P only in very low-density 
regions (where the rest-mass density po < lOOpatm) or 
when i/j^ > V^max- We set V^^ax = most simula- 

tions, except Case E where we need to set V^^ax = ^0 
stabilize the evolution. We note that even for this rela- 
tively small V^max^ ^ V^max holds truc only deep inside 
the AH, so this readjustment should not affect the dy- 
namics outside the BH. By relaxing this artificial pressure 
ceiling, we now observe a small remnant disk in many 
cases and achieve significantly better angular momentum 
conservation. For example, for Case A we previously had 
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5 J ^ 13% but now 6 J ^ 2% after relaxing the pressure 
ceiling. 



E. Effects of black-hole spin 

In Cases A-C, we vary the initial spin parameter of 
the black hole (parallel to the orbital angular momentum) 
from Jbi{/M^^=—0.5 to 0.75, fixing q = 3 and the initial 
orbital angular velocity MQ ^ 0.033. Notice in Table [J 
that as the initial BH spin parameter d increases, the 
total initial angular momentum increases, requiring more 
gravitational- wave cycles to emit angular momentum and 
bring the BH and NS close enough to merge. Thus we 
expect the binary to undergo more orbits before merger 
as Jbh/^bh increases. This is precisely what we find: for 
spins —0.5, 0.0, and 0.75, the binary inspiral phase lasts 
for 3.25, 4.5, and 6.5 orbits, respectively. Figure [9] shows 
the coordinate trajectories of the AH and NS centroids 
for Cases A and B. 

Figure \W\ shows snapshots of rest-mass density con- 
tours for Case A. The upper panels show the binary at 
the start of the simulation, two orbits into the simula- 
tion, and 3.5 orbits into the simulation, corresponding 
to the time at which the first few density contours have 
crossed into the AH. Notice that the equilibrium shape 
of the NS is for the most part undisturbed after 2 orbits. 
The lower panels in Fig. [10] are snapshots taken at ^44, 
^74, and ^96- They demonstrate how the NS tail deforms 
into a quasistationary disk, as the bulk of the matter is 
accreted onto the hole. 

Similarly, the upper plots of Fig. [TT] demonstrate that 
for Case B, the initial NS (upper left) retains its shape 
after 4 orbits (upper middle), and begins shedding its 
outer layers due to tidal disruption at about 5 orbits 
(upper right). Notice however, in this case, the tail at 
teo (lower left) is quite massive, and at ^74 (lower center) 
it is much larger than the tail at ^74 in Case A (lower left 
plot in Fig. [To]). The lower right plot is a snapshot taken 
near the end of the simulation, when a quasistationary 
disk resides outside the BH, at t = ^85 = 1878. 5M. 

To compare the effects of BH spin for all cases in this 
study, we plot the density contours at ti and tso for 
Cases A-C in Fig. [121 In the d = —0.5 case (Case C), 
the NS is basically "swallowed whole" by the BH during 
merger, leaving < 1% of the NS rest mass as a disk. As 
the spin increases, NS tidal disruption becomes more pro- 
nounced, resulting in long tidal tails that eventually form 
disks with rest mass ^ 4% and ~ 15% the rest mass of 
the NS in Cases A and B, respectively (see Fig. [13]). This 
result is not surprising. The ISCO decreases as the BH 
spin increases, and hence the tidal disruption of the NS 
occurs farther from the ISCO as the BH spin increases. 
Thus BHs with higher spin would likely lead to a more 
massive disk. 

Figure [13] shows the rest mass outside the AH as a 
function of time. We see that there are two phases of 
matter falling into the BH: a plunge phase and an accre- 



tion phase. The plunge phase occurs early in the merger 
as part of the NS matter streams onto the BH and the rest 
deforms into a tidal tail. The plunge phase ends when 
70%-90% of the NS matter falls into the BH, resulting in 
a sudden drop in the slope of the exterior Mq vs. time plot 
in Fig. [13] The matter in the tail, having larger specific 
angular momentum, spreads out and forms a disk. Ma- 
terial with lower specific angular momentum in the disk 
accretes onto the BH. Since there is neither viscosity nor 
magnetic fields in our simulation, the accretion should 
eventually cease as the evolution continues. However, in 
realistic astrophysical environments, magnetic fields do 
exist and could have substantial infiuence on the dynam- 
ics of the disk on secular timescales. In addition, some 
material in the disk is shock heated to high temperature 
(see Sec. IV Gp . Copious amounts of neutrinos and an- 
tineutrinos may be generated as a result. Both magnetic 
fields and neutrino-antineutrino pairs could drive ener- 
getic jets, resulting in an SGRB (see Sec. [VG] for further 
discussion). 



F. Effects of varying the mass ratio 

Figure [14] demonstrates how different mass ratios alter 
the NS material falling into the BH. For the q = 1 and 
q = 3 cases (Cases E and A), we clearly identify the 
plunge and accretion phases mentioned in the previous 
subsection. However, for the q = 5 case (Case D), the 
NS basically plunges into the BH, leaving < 1% of its 
rest mass outside the BH at the end of the simulation. 
This result is not surprising since, at the ISCO, the tidal 
effect of the BH is smaller for larger resulting in tidal 
disruption occurring closer to the ISCO. Moreover, for a 
fixed NS compaction, the horizon size of the BH is larger 
for larger q. Hence more NS matter is expected to fall 
into the BH. 

Figure [15] shows snapshots of density and velocity pro- 
files of the three cases at ti and tso- We again see differ- 
ent structures of the NS tail among the three cases. The 
merger of Case E (g = 1 case) is particularly interesting. 
Figure [16] shows snapshots of rest-mass density and K 
contours. Since the BH is less massive in Case E, tidal 
disruption occurs at a farther binary separation (upper 
middle plots). The disrupted NS matter curls around the 
BH, forming a hot, low-density spiral (upper right plots) 
that winds around the AH and smashes into the tidal tail, 
generating a large amount of shock heating. Some of the 
heated NS matter transfers angular momentum to the 
other part of the tail and falls into the BH. The remain- 
ing matter in the tail deforms into an inhomogeneous disk 
(lower center plots) before settling into a quasistationary 
state, in which a high density, relatively low temperature 
torus of matter surrounds the BH (lower right plots). 
Analysis of the accretion versus time (Fig. [T4|) demon- 
strates that the plunge phase is relatively long compared 
to Case A. Because of the complicated interaction of the 
disrupted NS matter in the merging process, the disk 
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FIG. 9: Trajectories of the BH and NS coordinate centroids for Cases A (left) and B (right). The BH coordinate centroid 
corresponds to the centroid of the AH, and the NS coordinate centroid is computed via: XI = (J^ x^p^d^x)/(J^ p^dPx)^ where 
V is the total simulation volume. Note that this equation contains a different V than Eq. (|18p . 
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FIG. 10: Snapshots of density and velocity profiles at selected times for Case A. The contours represent the density in the 
orbital plane, plotted according to po = Po,maxl0~°'^^"'~° °^, (i=0, 1, ... 12), with darker greyscaling for higher density. The 
maximum initial NS density is Kpo,ma^ = 0.126, or po,max = 9 x lO^^g cm~^(1.4M0/Mo)^. Arrows represent the velocity field 
in the orbital plane. The black hole AH interior is marked by a filled black circle. In cgs units, the initial ADM mass for this 
case is M = 2.5 x lO"^(Mo/1.4M0) s= 7.6(Mo/1.4M0)km. 



mass is less than that of Case A. 



G. Disk profile 



In Cases A, B, and E, a substantial disk forms after the 
merger. Figure [TTl shows the density and thermal energy 
profiles of the disk in both the equatorial and meridional 




FIG. 11: Snapshots of density and velocity profiles at selected times for Case B. The contours represent the density in the 
orbital plane, plotted according to po = PcmaxlO"^"^^-^"^ ^'^ 0=0, 1, ... 12), with darker greyscaling for higher density. The 
maximum initial NS density = 0.126, or /9o,max — 9 X 10"'^^g cm ^(1.471^0/71^0)^- Arrows represent the velocity field 

in the orbital plane. We specify the black hole AH interior in each snapshot with a filled black circle. In cgs units, the initial 
ADM mass for this case is M = 2.5 x lO"^(Mo/1.4M0)s= 7.6(Mo/1.4M0)km. 



plane at the end of the simulation for Cases A, B, and E. 
For ease of comparison, we specify the length scale in km, 
assuming the NS has a rest mass of I.4M0; it can be con- 
verted to units of M via the formula M = 1.9{q -\- l)km. 
The most distinguishable features of the disks are their 
characteristic radius raisk and final mass, Mdisk- The 
characteristic disk radii lie between 20km and 100km, 
as seen in the top row of plots. The q = 1 mass ratio 
(Case E) produces the smallest (rdisk = 20km), densest 
(maximum density po,max ~ 6 x lO^^g cm~^), and least 
massive disk (Mdisk/^o ^ 2%). Meanwhile, the result- 
ing disk in the canonical q = 3 case (Case A) is about 
2.5 times larger than Case E (rdisk = 50km), but pos- 
sesses a lower maximum density by an order of magnitude 
(?^ 4 X lO^^g cm~^). Due in part to its larger volume, the 
mass of the disk in Case A is about 70% larger than in 
Case E (Mdisk/^o ^ 4%). The maximum density of 
Case B's disk is po,max ~ 5 x lO^^g cm~^, which is simi- 
lar to Case A. However, the disk is about twice as large 
in size, with a characteristic radius of rdisk = 100km, 
yielding a total disk mass of Mdisk /^o ^15%. 

Despite their differences, the disks in Fig. [T71 share 
many common features. For example, the top row of 
plots show that the disk forms a torus whose density 



plummets at the BH ISCO (roughly at coordinate radius 
7km, 25km, and 12km for Cases E, A, and B, respec- 
tively, in the equatorial plane H^). Also, the characteris- 
tic height of the disks is about 15%-20% of the character- 
istic radius in each case (top two rows of plots). Finally, 
the higher density regions tend to be colder than lower 
density regions, but K ^ 85 everywhere in the disks in 
Cases A and B, and K ^ 6 everywhere in Case E's disk. 
The characteristic K in the disks in Cases A, B, and E 
are roughly 200, 200, and 10, respectively. 

As in Paper I, the temperature T can be estimated 
roughly by modeling the temperature dependence of the 
specific thermal internal energy density eth as 



3kT 
2mn 



f 



Po 



(20) 



(compare [83), where rUn is the mass of a nucleon, k is 
the Boltzmann constant, and a is the radiation constant. 
The first term represents the approximate thermal en- 
ergy of the nucleons, and the second term accounts for the 
thermal energy due to radiation and (thermal) relativistic 
particles. The factor / reflects the number of species of 
ultrarelativistic particles that contribute to thermal radi- 
ation. When T <C 2me/k ~ lO^^K, where rrie is the mass 
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FIG. 12: Snapshots of density and velocity profiles at ti (top row) and tso (bottom row) for Cases C, A, and B, which have 
initial BH spins of a = —0.5, 0.0, and 0.75, respectively. The contours represent the density in the orbital plane, plotted 
according to po = Po,maxl0~°"^^"'~°"°'^ (^=0, 1, ••• 12), with darker greyscaling for higher density. The maximum initial NS 
density is /^po,max = 0.126 for all cases, or po,max = 9 x lO^^g cm~^(1.4M0/Mo)^. Arrows represent the velocity field in the 
orbital plane. The black hole AH interior is marked by a filled black circle. In cgs units, the initial ADM mass for these cases 
is M = 2.5 X lO"^(Mo/1.4M0)s= 7.6(Mo/1.4M0)km. 



of electron, thermal radiation is dominated by photons 
and / = 1. When T ^ 2me//c, electrons and positrons 
become ultrarelativistic and also contribute to radiation, 
and / = 1 + 2 X (7/8) = 11/4. At sufficiently high tem- 
peratures and densities (T ^ lO^^K, po ^ 10^^ g cm~^), 
thermal neutrinos are copiously generated and become 
trapped, so, taking into account three flavors of neutri- 
nos and anti- neutrinos, / = 11/4 + 6 x (7/8) = 8. 



The characteristic density and K in the disk for 
Cases A and B are npQ ~ 4 x 10~^ (or po ~ 3 x 
IQiig cm-3) and K - 200. We flnd T - 5 x lO^^i^ 
(or ~4MeV). For Case E, we have hipQ ~ 2 x 10~^ (or 
po lO^^g cm-3). We flnd T - lO^^K (or -IMeV) for 
Case E's disk. According to numerical models of rotating 
BHs with ambient disks in [43^, the remnant disks could 
produce a total gamma-ray energy E ~ lO^^-lO^^erg 
from neutrino-antineutrino annihilation. This result is 
promising for generating a SGRB. 



H. Gravitational wave emission 

Following the literature, we decompose the gravita- 
tional waveform /i+ and into s = —2 spin- weighted 
spherical harmonics -2X^771 as follows 

h+ - ih^ = - ih'T) -2Yim , (21) 

1,171 

where h^J^ and are real functions. Each (/,m) mode 
is a function of radius and time only. We extract grav- 
itational waves using both the Regge- Wheeler- Zerilli- 
Moncrief formalism and the Newman-Penrose Weyl 
scalar 1^4. We flnd good agreement between waveforms 
computed by both methods, as demonstrated in Paper I. 
Here we only present the waveforms computed from i/j^. 

Gravitational waves are extracted at several different 
radii. We flnd that the extracted waveforms overlap ex- 
tremely well provided the extraction radius rex ^ 40M, 
as shown in Fig. [T8l 

Figures [T9l and [2Ql demonstrate that the maximum am- 
plitude of the waves is similar in Cases A, B, C, D, and 
E. With approximately the same initial MQ^ we observe 
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FIG. 13: Rest-mass fraction outside the BH for different ini- 
tial BH spins (Cases C, A, and B). Here, the time coordinate 
is shifted by t25, the time at which 25% of the NS rest mass 
has fahen into the apparent horizon. 



more wave cycles as the BH spin is increased (Fig. [T9|) 
and as the mass ratio q is decreased (Fig. [20]), consistent 
with the discussion in Sees. IV El and IV Fl Table [TTll spec- 
ifies the radiated energy AEqw and angular momentum 
A Jgw in each case, as well as the kick velocity 'Ukick due 
to recoil. 



For Case E, NS disruption occurs at relatively low or- 
bital frequency. In fact, the total radiated angular mo- 
mentum AJgw reaches 90% of its final value at tio (the 
time at which only 10% of the NS has been accreted), 
and 99% of its final value at tso- Hence there is a large 
amount of angular momentum still remaining in the NS 
matter before it is accreted, and during the accretion 
phase about 98% of the NS rest mass is accreted. Thus 
the bulk of the orbital angular momentum is used to 
spin up the BH. The final disk mass is only ~ 1% of 
the binary's ADM mass, so we can reliably estimate the 
final value of the BH spin parameter using the Kerr for- 
mula (Eq. (p!6|)). We obtain a = 0.85, which is signifi- 
cantly larger than all other cases that result in a simi- 
larly small disk. Compared to the final spin parameter 
from equal-mass BHBH mergers, a = 0.686 (e.g. [89] and 
Table [VH] in Appendix [X]) , we find that matter accretion 
in BHNS binary mergers is much more effective at spin- 
ning up black holes. This is not unexpected, since NS 
tidal disruption spreads the NS matter around the BH, 
diminishing the multipole moments of the system, and 
the associated gravitational wave emission. 
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FIG. 14: Rest-mass fraction outside the BH for different 
BHNS mass ratios. Here, the time coordinate is shifted by 
t25 , the time at which 25% of the NS rest mass has fallen into 
the apparent horizon. 



Figures [2T1 and [221 show the effective gravitational wave 
strains in frequency space for all cases, computed via the 
method summarized in Sec. HI F of Paper I, and in- 
cluding only the dominant (2, 2) and (2, —2) modes. For 
comparison against the Advanced LIGO sensitivity curve 
hLiGoif) = \/fSUf) ([90]), we plot BHNS wave strains 
for the cases we study, assuming the NS has a rest mass of 
IAMq and binary distance of lOOMpc. This is the dis- 
tance required to reach one merger per year, assuming 
an overall rate of 10 mergers per Myr per Milky V^y- 
equivalent galaxy (and a density of 0.1 gal/Mpc^) jgij . 
This distance is roughly that of the Coma cluster, and ap- 
proximately five times the distance to the Virgo cluster. 
The gravitational wave spectra of nonspinning BHBH 
mergers (computed from Eqs. (4.12)-(4.19) of [88]) with 
the same mass ratio as in Cases E, A and D are plotted 
in the same graphs for comparison. We see that the wave 
signal drops significantly at and above the frequency cor- 
responding to the onset of NS tidal disruption. The dif- 
ference in wave signals between BHBH and BHNS merg- 
ers is marginally observable in the Advanced LIGO fre- 
quency band in most cases. Distinguishing BHNS from 
BHBH inspirals and mergers may require narrow-band 
detection techniques with advanced detectors. The ob- 
servation of an accompanying SGRB would also serve to 
distinguish the events. 
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FIG. 15: Snapshots of density and velocity profiles at ti (top row) and tso (bottom row) for Cases E, A, and D, which 
have initial mass ratios of ^ = 1, 3, and 5, respectively. The contours represent the density in the orbital plane, plotted 
according to po = po,maxlO~°"^^"'~°"°'^ (^=0, 1, ••• 12), with darker greyscaling for higher density. The maximum initial NS 
density is Kpo,niax = 0.126 for all cases, or po,max = 9 x lO^^g cm~^(1.4M0/Mo)^. Arrows represent the velocity field in the 
orbital plane. The black hole AH is marked by a filled black circle. In cgs units, the total ADM mass for these cases is 
M = 6.25 X 10"^(^ + l)(Mo/1.4M0)s= 1.9{q + l)(Mo/1.4M0)km. 



VI. SUMMARY AND DISCUSSION 

In this paper we perform a new set of fully general rela- 
tivistic BHNS simulations using our code upgraded with 
AMR capability. We vary the initial binary separations 
and confirm that the simulation outcomes do not change 
when the initial BHNS separation is large enough. In 
most cases, a toroidal disk forms around the black hole 
remnant after the merger. This result revises our find- 
ing of Paper I, where we artificially imposed a uniform 
pressure ceiling on all the matter present in the simula- 
tion. This ceiling removed a substantial amount of angu- 
lar momentum in the NS matter during the merger phase, 
causing the matter to fall into the BH prematurely, and 
suppressing disk formation. 

Fixing the mass ratio (7 = 3, we study the effects of BH 
spin by evolving initial data with BH's spin parameter 
a = —0.5 (count-rotating), and 0.75 with nearly the 
same initial orbital angular velocity MQ ^ 0.033. Not 
surprisingly, we find that the BHNS inspiral phase lasts 
longer as d increases, and the final disk mass increases 
from < 1% (a = -0.5) to ^ 4% (a = 0) and ^ 15% 



(a = 0.75) of the NS's rest mass. 

Next, we study the effect of varying the mass ratio q 
by evolving BHNS initial data with q = 1 and 5. For 
the q = 5 case almost all the NS matter plunges into 
the BH, leaving < 1% of the NS's rest mass to form a 
disk at the end of the simulation. For the q = 1 case, 
a low-density, hot spiral region of disrupted NS matter 
winds around the BH and smashes into the nascent NS 
tidal tail, causing strong shock heating before it rapidly 
falls into the BH. Most of the remaining matter in the 
tail eventually falls into the BH, leaving about 2% of the 
NS's rest mass behind to form a disk. 

The disks formed after the merger of Cases A and B 
have a characteristic density of 3 x lO^^g cm~^ and a 
temperature of T ~ 5 x lO^^K, assuming the NS's rest 
mass is Mq = 1.4M0. A pply ing the results of the sim- 
ulations of BH disks in [43], the disk may produce a 
gamma-ray energy of ^ ^ 10^^-lO^^erg due to neutrino- 
antineutrino annihilation. Hence the merger of a BHNS 
binary is a promising central engine for a SGRB. But self- 
consistent simulations, taking into account the detailed 
microphysics and including the correct EOS, are neces- 
sary to fully explore the viability of and the mechanism 
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FIG. 16: Snapshots of rest-mass density po and K contours for Case E = 1 mass ratio). The first and third rows show 
snapshots of density contours and velocity profiles in the orbital plane. The second and fourth rows show snapshots of 
K — P/kpq. The density is plotted according to po = yOo,maxlO~°'^^'^~°"°^ 0=0, 1, ••• 12), with darker greyscaling for higher 
density. The maximum density of the initial NS IS /Tj/?0,max = 0.126, or po,max = 9 X 10"^^g cm ^{IAMq/Mq)'^ . Arrows represent 
the velocity field in the orbital plane. K is plotted according to K = ]^Q0-32j (^j—q^ 12), with darker greyscaling for higher 

K. These figures neglect K in regions where the density is less than po,min. The black hole AH interior is marked by a filled 
black circle. In cgs units, M = 1.3 x lO"^(Mo/1.4M0)s=3.8(Mo/1.4M0')km. 
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FIG. 17: Snapshots of rest-mass density po^ K — P/kpq, and velocity profiles at the end of the simulation, for Cases E 
(left column), A (middle column), and B (right column) on the orbital (xy, first and third rows) and meridional {xz and yz, 
second and fourth rows) planes. In the top six plots, density contours are shown according to po = po,maxlO~°"^-^~^"^ (i=0, 
1, ... 13) for Cases A and B, po = po,maxlO~°-^^-^~^-^^ (i=0, 1, ... 12) for Case E. The maximum density of the initial NS is 
= 0.126 for all cases, or po,max — 9 x 10''^^g cm "^(1.4Ai^0/Aio)^- The bottom six plots show contours of plotted 
according to K = 10°-^^^' (j=0, 1, ... 12). The ii^ fi gures neglect regions where the density is less than po,min. All contour plots 
show darker grey scaling for higher po and K. Arrows represent the velocity field in the given plane. The AH interior in the 
orbital plane is marked by a filled black circle. Length scales are specified in km, assuming the NS has a rest mass of IAMq; 
it can be converted to units of M via the formula M = 1.9(^ + l)km. 
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FIG. 18: Gravitation wave signal from Case A. Shown here 
are the (2,2) and (3,3) modes of r/i+(t — r) and rhx{t — r) 
extracted at radii rex = 43AM (black solid lines), 54. 2M (red 
dotted lines) and 83. 2M (blue dashed lines). Note that the 
three lines almost overlap in each case. 
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FIG. 19: Gravitational wave signal from (top to bottom) 
Cases C, A and B. Black solid (blue dash) line denotes the 
(2,2) mode of /i+ {h^) extracted at rex = 43AM. 



FIG. 20: Same as Fig. [19] but for cases (top to bottom) E, A 
and D. 



for generating SGRBs from BHNS mergers. 

We compute the BHNS gravitational waveforms and 
power spectra, and find that the signals are attenuated 
at frequencies roughly equal to double the orbital fre- 
quency at which tidal disruption begins. The resulting 
deviation of BHNS and BHBH waveforms is visible in the 
Advanced LIGO high frequency band out to distances 
^ 100 Mpc. Should the chirp mass determination, com- 
bined with higher order PN waveform phase effects, al- 
low for an independent determination of the component 
masses and spins of the binary companions, the measure- 
ment of the BHNS tidal disruption frequency should give 
a good estimate of the NS radius and, hence, insight into 
the nuclear EOS. 

Carpet-based AMR has provided us with a great im- 
provement in efficiency and accuracy, and future work 
will continue in this direction. To improve code perfor- 
mance (which is usually memory-bound), we also plan 
to employ new techniques (e.g., loop-tiling) for min- 
imizing cache misses. To improve accuracy, we will 
implement the piecewise parabolic method (PPM) or 
weighted essentially non-oscillatory (WENO) reconstruc- 
tion schemes in hydrodynamics/MHD and experiment 
with tapered grids to replace second order temporal pro- 
longation (see 0). We plan to study techniques for min- 
imizing refinement boundary crossings, by adding more 
refinement boxes around the disrupted NS. This might 
help reduce our angular momentum conservation viola- 
tion to below 1% for all cases. 

Rich with the challenges of modeling extreme matter 
in the presence of intense gravitational fields, BHNS bi- 
nary mergers will likely remain a topic at the forefront 
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FIG. 21: Gravitational wave spectrum for the Cases C, A and 
B, computed via the method summarized in Sec. IIIF of Pa- 
per I. The sohd curve shows the combined waveform found 
by attaching the restricted 2.5 order PN waveform to our nu- 
merical signal (including only the dominant (2, 2) and (2,-2) 
modes), while the dotted curve shows the contribution from 
the latter only. The dashed curve is the analytic fit derived by 
[ssj from analysis of multi-orbit nonspinning BHBH binaries 
with the same mass ratios q as the BHNS, which maintain 
significantly more power at higher frequencies. The heavy 
solid curve is the effective strain of the Advanced LIGO de- 
tector, defined such that /iligo(/) = yf^hXf)- To set phys- 
ical units, we assume a NS rest mass of Mo = IAMq and a 
source distance of i^=100Mpc. 



of theoretical astrophysics and numerical relativity for 
years to come. The most exciting possibility remains the 
simultaneous detection of a gravitational wave signal and 
a ORB. Our theoretical work on BHNSs is partly in an- 
ticipation and preparation for that discovery. 
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FIG. 22: Same as Fig.[2T]but for Cases E, A and D. 



APPENDIX A: DETAILS OF BHBH TEST 
SUMULATIONS 

Simulating the late inspiral, merger and ringdown of 
a binary black hole was for many years the "holy grail" 
of numerical relativity. Dozens of researchers have spent 
many years attempting to formulate a stable algorithm 
capable of solving this important problem. In this ap- 
pendix we summarize two BHBH simulations in detail. 
Performing a simulation of binary black hole coalescence 
provides a highly nontrivial and reasonably comprehen- 
sive laboratory for testing 3+1 codes that solve Einstein's 
vacuum field equations in the strong-field regime of gen- 
eral relativity. It tests the caliber of not only the ba- 
sic evolution scheme, but also the all-important diagnos- 
tic routines. These routines measure globally conserved 
quantities like the mass and angular momentum of the 
system, the location and measure of all black hole hori- 
zons, the asymptotic gravitational waveforms, the recoil 
motion of the black hole remnant, etc. Appreciable ef- 
fort must be expended to implement reliable diagnostic 
routines in order to extract useful physical information 
from the numerical output. 

We consider two cases of merging, nonspinning black 
holes to test our code: an equal- mass {q = 1) and an 
unequal-mass {q = 3) system. Both cases are constructed 
via the binary puncture technique discussed in [74] , with 
parameters as listed in Table HVl Our code employs AMR 
with equatorial symmetry using the moving-box Carpet 
infrastructure, so that the boxes track the AH centroids 
throughout the simulations. The grid setup is described 
in Table [Vl The basic evolution scheme and numerical 
techniques are described in Sees. [Ill and Hill A review 



23 



of the special algorithms and parameters chosen in these 
simulations is summarized in Table IVIi and the results 
are summarized in Table IVlTl The coordinate trajectories 
of the two black holes are plotted for the two cases in 
Fig. [23l and the gravitational waveforms are plotted in 
Fig. El 

The results found here are in good agreement with 
those reported by previous investigators [50, 77, 8o!, [8l|. 
Of special significance is the degree to which total energy 
and angular momentum are conserved when properly ac- 
counting for losses due to gravitational wave emission. 
The fractional errors are seen in Table IVIII to be a few 
times 10~^ for the energy and 10~^ for the angular mo- 
mentum. Also, different measures of the mass and spin 
of the final (stationary) remnant Kerr black hole, such as 
the irreducible mass and the ADM mass, agree closely, 
given the adopted computational resources. 



APPENDIX B: SHOCK HEATING IN 
LOW-DENSITY REGIONS 

In this appendix, we estimate the amount of thermal 
energy generated by strong shocks in low-density regions. 
Such shocks are generated when fluid elements of initially 
cold matter collide. For simplicity, we assume Newtonian 
hydrodynamics in our analysis, which is adequate for a 
rough estimate, since the fluid velocity du ring t he BHNS 
merger is only mildly-relativistic {v ~ ^yM/r ^0.4 for 



> 



6M) and P <C po in low-density regions. 



Let pi. Pi, vi, and hzi = Pi/Pi be the upstream rest- 
mass density, pressure, velocity, and poly tropic constant, 
respectively as measured in the shock frame, and p2^ P2, 
V2 and K2 the corresponding quantities in the downstream 
fluid. The upstream Mach number is given by 



M 



Vi 
Cl 



rPi 

pi 



-1/2 



Vl 



r-i 



(Bl) 



where ci is the sound speed in the upstream fluid. It can 
be shown that (see e.g. Sec. 89 of [94]) 



P2 

Pi 
P2 
Pi 



(r-i)A^2 + 2 ' 

2TM'^ r - 1 



(B2) 
(B3) 



In the shock frame, the flow always enters the front su- 
personically, (i.e. M > l)^ hence it follows from Eqs. ()B2p 



and ()B3p that P2/P1 > 1 and P2/P1 > 1. Also, for adia- 
batic flow on either side of the shock, 



Pi 

where K' 
obtain 



K' = 



Pi 
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^2/^1- Combining Eqs. ()B3p and ()B4p . we 



2rA^2 - (F - 1) 



1 



2 



(B5) 



(r + l)r+i 

As a shock generates entropy, we must have K' > 1. 
In fact, it is easy to prove from Eq. ()B5p that K' > 1 
whenever M. > 1 and F > 1. The thermal energy of the 
downstream fluid can be measured by = P2I = 
(k\Ik)K' ^ where n is the polytropic constant for a cold 
gas. Since we always have k\ > k (with the equality 
holding for cold, unshocked matter), we must always have 
K > 1 when Ai > 1. In other words K exceeds unity 
across a shock front. Using Eqs. (|Bip . (|B5p and K = 
(Ai:i/Ai:)i^', we have 



K 



2vl 



(F + \)kp\ 



r-i 



(B6) 



in the strong shock limit (A^ ^ 1). Equation (|B6P 
shows that when the density p\ of the upstream mat- 
ter is low, the amount of the shock heating downstream 
can be substantial, i.e. K ^ \. The initial maximum 
rest-mass density of the NSs in all of our BHNS models 
is /^/)o,max = 0.126. Setting F = 2, we have 



K , 



V0.4 



Pi 



(B7) 



in the strong shock limit. We see that imposing the pres- 
sure ceiling P < lO/^pg (i-^- ^ ^ ^0) spuriously 
remove a substantial amount of thermal energy and pres- 
sure generated by strong shocks when the density is below 
0.01po,max- Such low-dcusity regions contribute a sub- 
stantial fraction of the NS rest mass during the merger 
phase following tidal disruption. The pressure ceiling 
then serves to remove a substantial amount of torque 
and angular momentum from the fluid, causing a large 
violation of angular momentum conservation. 
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TABLE IV: Initial data of our BHBH simulations 





Equal mass 


Unequal mass 


Location of punctures 


£c+ = (0,4.891,0), 


£c+ = (5.25,0,0), 




x_ = (0,-4.891,0) 


X- = (-1.75,0,0) 


"Bare" mass 


m+ = m_ = 0.4856 


m+ = 0.234, m- = 0.735 


Spin 






Momentum 


P+ = (0.0969,0,0), 


P+ = (0,0.09407, 0), 




P- (-0.0969, 0,0) 


p_ = (0, -0.09407, 0) 


ADM mass of the system 


M = 0.9894 


M = 0.9895 


ADM angular momentum 


J = 0.9479 


J = 0.6587 


Irreducible mass of the BHs 


M+ = M-= 0.5000 


M+ = 0.2498, M- = 0.7494 


Mean coordinate radius of BHs' horizon 


7^+ = 7^- = 0.2359 


7^+ = 0.1097, 7^_ = 0.3605 


Binary coordinate separation ( x+ — X-\) 


9.887M 


7.074M 



TABLE V: AMR grid setup in our BHBH simulations. 





equal- mass 


unequal mass 


L 


320 (323.4M) 


409.6 (413.9M) 


N+, N- 


8, 8 


10, 9 


R+ 


(160, 80, 140/3, 20, 5, 2.5, 1.25, 0.625) 


(256, 128, 64, 24, 12, 6, 3, 1.5, 0.75, 0.375) 


A+ 


(4, 2, 1, 0.5, 0.25, 0.125, 0.0625, 0.03125) 


(5.12, 2.56, 1.28, 0.64, 0.32, 0.16, 0.08, 0.04, 0.02, 0.01) 


R- 


(160, 80, 140/3, 20, 5, 2.5, 1.25, 0.625) 


(256, 128, 64, 32, 16, 8, 4, 2, 1) 


A- 


(4, 2, 1, 0.5, 0.25, 0.125, 0.0625, 0.03125) 


(5.12, 2.56, 1.28, 0.64, 0.32, 0.16, 0.08, 0.04, 0.02) 




^ 15,15 


^ 22,36 



L: Location of the outer boundary Xmax — 2/max - Zma^ - L, Xmin - 2/min - -L, Zmin - 0. 

iV+ (N-): Number of refinement levels centered at the "+" ("-") puncture. 

(R-): Half side length of refinement boxes centered at the ("-") puncture (in code unit). 
A+ (A-): Grid spacing in each refinement level centered at the ("-") puncture (in code unit). Note that the grid spacing 
doubles in each successive refinement level. The coarsest grid spacing is twice of that of the outermost refinement level, which 
is 8 for the equal-mass case and 10.24 for the unequal-mass case. 

^AH (^ah)- Number of grid points across the mean diameter of the apparent horizon of the ("-") BH at t = 0. 




FIG. 23: Trajectories of the coordinate centroid of the "+" (black solid line) and "-" (blue dashed line) black hole from the 
equal-mass (left panel) and unequal-mass (right panel) BHBH simulations. 
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TABLE VI: Evolution technique used in our BHBH simulations 



Formalism 


BSSN-based 


Lapse 


1+log slicing: = —2aK, do = dt — P^dj. 


Shift 


"Gamma-freezing": dof3' = (3/4)B\ doB' = - r]B\ r] = 0.25 = 0.2474/M. 


Conformal variable 


Evolve W = e~^'^ instead of 0. 


Symmetry imposed 


Equatorial (i.e., symmetry across equatorial plane). 


Other Numerical techniques 


(1) Enforce auxiliary constraints det(7ij) = 1 and Ti (Aij) = using the technique 
described in Paper L 

(2) Enforce auxiliary constraint F* = — ^j7*"^ using the technique as in, e.g. [93]. 

(3) Add 5th order Kreiss-Oliger dissipation of the form 
ie/64){Ax^d^^ + Ay^dl + Az^dt)f 

for all BSSN, lapse and shift variables /. Here Ax, Ay and Az are grid spacings 
and we set the strength e = 0.1. 


Grid-driver code 


Carpet 


Spatial differencing 


4th order upwinded on shift advection terms, 4th order centered everywhere else 


GW extraction radii 


10.11M-54.58M for the equal-mass case, 30.32M-70.75M for the unequal-mass case. 


Temporal differencing 


MoL: 4th order Runge-Kutta (RK4) 


Courant-Friedrichs-Lewy (CFL) factor 


0.25 


Prolongation 


5th order spatial prolongation, 2nd order temporal prolongation 
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TABLE VII: Summary of our BHBH simulation results 







nnpnnpil m?i<i;<i; ( n — ^ i 


7\ /T / 7\ /f 

Mbh/M 


u.yoi ( 


n nonn 

u.youy 


Jbh/Mbh 


U.Dooz 


U.o41o 


AEgw/M 


0.03794 


0.01934 


AJgw/M^ 


0.3306 


0.1513 


'i^kick 




174 km s"^ 




4 X 10"^ 


-2 X 10"^ 


5J 


4 X 10"^ 


9 X 10"^ 



Mbh and Jbh/M^hI Mass and spin parameter of the final 
BH, as determined by its irreducible mass and the ratio of 
polar and equatorial circumferences (using the Kerr formula, 
given in Eq. (|16p V Here M is the initial ADM mass as listed 
in Table HVl 

AEgw and AJgw: Energy and angular momentum carried 
off by gravitational radiation. GWs are extracted at 50.53M 
using '04- 

'^^kick^ The kick velocity measured by GW. Following 
Gonzalez et al jSll], GW data with t-r < 50M, 
corresponding to the "junk" GW, for the unequal-mass case 
are removed before computing the kick velocity. 
SE = (M — Mbh — AEqw) / M is a measure of violation of 
energy conservation. 

5J = (J — Jbh — AJgw) /J is a measure of violation of 
angular momentum conservation. 








300 600 



900 



C\2 




100 200 300 

(t-r)/M 

FIG. 24: Gravitational waveform from the equal-mass (up- 
per panel) and unequal-mass (lower panel) BHBH merger. 
Shown here are {r/M)hf (black solid line) and {r/M)hf 
(blue dashed line) as a function of the retarded time t — r 
extracted at radius r = 50.53M. 
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